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Singularly perturbed phase response curves 
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Abstract —In this paper we propose a novel geometric 
method, based on singular perturhations, to approximate 
isochrones of relaxation oscillators and predict the qualitative 
shape of their (finite) phase response curve. This approach 
complements the infinitesimal phase response curve approach 
to relaxation oscillators and overcomes its limitations near the 
singular limit. We illustrate the power of the methodology hy 
deriving semi-analytic formula for the (finite) phase response 
curve of generic pianar retaxation osciltators to impulses 
and square pulses of finite duration and verify its goodness 
numerically on the FitzHugh-Nagumo model. 

I. Introduction 

The phase response curve (PRC) characterizes the input- 
output behavior of oscillatory systems [1], [2]. It has wide 
applications ranging from oscillator control [3], [4] to the 
analysis of oscillator network synchronization [5], [ 6 ]. Sys¬ 
tematic and analytic prediction of an oscillator phase response 
curve is a hard task in general and it can be accomplished only 
in very specific cases. This usually leads to intense numerical 
investigations, which might weaken the relevance of phase 
response curve approach in engineering applications. 

The classical approach relies on analytically computing the 
infinitesimal (linearized) phase response curve and then use 
convolution to compute the phase response curve for generic 
inputs [7]. However, because the measure of the region where 
the linear approximation is valid shrinks to zero as the time- 
scale separation is increased, this approach is valid only for 
inputs of infinitesimally small amplitude [ 8 ]. 

To overcome this limitations, we used a complementary 
approach and study the global structure of the oscillator 
isochrones from a geometric viewpoint, using singular per¬ 
turbation method [9]. Based on this analysis we derived semi- 
analytic formulas to approximate the (finite) phase response 
curve to arbitrary inputs in the singular limit. As opposed 
to the infinitesimal phase response curve approach, the error 
between the real and the predicted singular phase response 
curve goes to zero as the time-scale separation increases. The 
method is thoroughly illustrated by predicting the phase re¬ 
sponse curve of a generic relaxation oscillator to impulses and 
square pulses of finite duration. These results are a first step 
toward a geometric theory for (finite) phase response curves 
of singularly perturbed oscillators, including bursters [ 10 ]. 

The results of the article are primarily drawn from the Ph.D. 
dissertation of the first author [ 11 ]. 
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This paper is organized as follows. Section [H] summarizes 
basics concepts of singular perturbation theory and describes 
the underlying geometry of relaxation oscillator dynamics. 
Sectionj^stresses the limitation of the infinitesimal phase re¬ 
sponse curve approach in the context of relaxation oscillators 
with inputs acting on the fast variable dynamics. Section m 
introduces the novel concept of singularly perturbed phase 
response curves predicted from the singular limit. Section [V] 
illustrates our geometric approach on the FitzHugh-Nagumo 
neural model. 

II. Relaxation oscillators and their geometry 

We consider a two-dimensional fast-slow dynamical system 
of the form 

x= f{x) — z + u, (la) 

z = eg{x,z), (lb) 

where ' denotes differentiation with respect to the time t, 
{x, z) € K^, u G M., and 0 < e ^ 1. The solution at 
time t to the initial value problem ([T]) from the initial condition 
{xo^zfi) G at time 0 is denoted by {xq, Zo),u{-)), 
with ^1(0, (xo, zq), u(-)) = {xo,zo). In the slow time scale 
T := et, dynamics Q become 

ex' = f{x) — z + u, ( 2 a) 

z'=g{x,z), (2b) 

where ' denotes differentiation with respect to the slow time r. 
For e ^ 0, 0 and 0 are equivalent. We call 0 the fast 
dynamics and 0 the slow dynamics. In the limit e —0, we 
obtain from 0 and 0 the layer dynamics 

X = f{x) — z + u, (3a) 

i = 0, (3b) 

describing the fast evolution far from the critical manifold 
:= {(t, z) G : f{x) — z + u = O}, and the reduced 
dynamics 

0 = f{x) — z + u, (4a) 

z'=g{x,z), (4b) 

describing the slow evolution along S'^. 

Under some mild technical assumptions [9, Theorem 2.1], 
in particular that the critical manifold 5° is S-shaped, the 
zero-input system 0 has a unique periodic orbit 7 *^ sliding 
along the stable branches of and shadowing the singular 
periodic orbit 7 ° illustrated in Figure The singular periodic 
orbit 7 ° is defined as the union of two pieces of the critical 
manifold associated with a slow evolution and two critical 
fibers associated with jumps. 
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Fig. 1. Geometry of relaxation oscillators. The critical manifold 
is a S-shaped curve. Under some technical assumptions [9], the singular 
system {TJ admits a singular periodic orbit 7 ® defined as the union of two 
pieces of the critical manifold associated with a slow evolution (green solid 
lines) and two critical fibers associated with jumps (green dashed lines). 

Remark 1: In the slow time scale, the singularly perturbed 
period T/ converges towards the singular period T°, which is 
equal to the finite time required to slide along both portions 
of the critical manifold (jumps are instantaneous), that is, 
lim£_j.oT!f =: T^. In the fast time scale, the singularly 
perturbed period Tf blows up to infinity, that is, limj_j.o =: 
T°, with lime_>.oI 7 = lime_>.o /e = +00. Note that 
corresponding angular frequencies are denoted uif = 27r jTf 
and Wj = 2Tr jT^. 

III. Limitation of the infinitesimal approximation 

FOR FINITE PHASE RESPONSE CURVES 

In this section, we introduce the concepts of phase map 
and phase response curves following the terminology of [1] 
and [2]. (The interested reader is referred to [7] for detailed.) 
We emphasize the limitation of the infinitesimal approxima¬ 
tion for finite phase response curves of relaxation oscillators. 

A. Phase map and isochrons 

Because of the periodic nature of its steady-state behavior, 
it is appealing to study the oscillator dynamics on the unit 
circle S^. The key ingredients of this phase reduction are the 
concepts of phase map and isochrons. 

The (asymptotic) phase map 0*^ : C —)■ is 

a mapping that associates with every point in the basin of 
attraction a phase on the unit circle S^. It is defined 

such that the phase variable 0(f) := (xq, 20 )) «(•)))> 

that is, the image of the flow through the phase map, linearly 
increases with time for the input 0. 

Isochron X'^(0) are the set of points in B{'-f‘^) that are 
associated with the same asymptotic phase 9 on the unit 
circle . Isochrons are level sets of the phase map. 

B. Phase response curves 

An input u(-) is phase-resetting if the solution of (0 forced 
by m(-) asymptotically converges to the periodic orbit. 
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Fig. 2. Qualitative trade-off between the infinitesimal approximation and 
the singular approximation as a function of the time-scale sepai'ation e. 

The (finite) phase response curve ^ 

[—TT, tt) associates with each phase the asymptotic phase shift 
of system 0 in response to a phase-resetting input u{-). 

The infinitesimal phase response curve —>■ K is the 

relative asymptotic phase shift of system 0 in response to an 
infinitesimal phase-resetting impulse (Dirac 5 function), that 
is, q\-) := lima^o 

C. Limitation of the infinitesimal approximation 

In the classical approach, the (finite) phase response curve 
is approximated by the “convolution” between the infinitesi¬ 
mal phase response curve and the phase-resetting input 

f+OO 

= / q^{uts + -)u[s)dsPO{\\u{-)\\^). 

Jo 

For relaxation oscillators, this approximation is only valid for 
inputs that are much smaller than the singular perturbation 
parameter, that is, 0 < ||m(-)II ^ e ^ 1 (see [8] for details). 
Therefore, the domain of validity of this approximation van¬ 
ishes in the singular limit (e —> 0). 

Intuitively, this limitation comes from the fact that, the 
singular trajectory induced by the input u{-) might jump 
instantaneously from one branch of the critical manifold to 
the other. This behavior involves a global phenomenon that 
cannot be captured by a local approximation. 

IV. Prediction from the singular limit 

FOR FINITE PHASE RESPONSE CURVES 

The main idea of our approach is to take advantage of time- 
scale separation to study the finite phase response curve in 
the singular limit. For a sufficiently small singular parameter 
0 < e ^ 1, the singularly perturbed finite phase response 
curve can be approximated by a singular finite 

phase response curve u{-)), that is, 

for any phase-resetting input u{-) and with 0 < ^ < 1. 
Geometric singular perturbation arguments let /3 ^ 1/2 [9]. 

Figure illustrates qualitatively the existing trade-off be¬ 
tween the infinitesimal approximation and the singular ap¬ 
proximation as a function of the time-scale separation e. 

In this section, we show how to exploit the geometry 
of relaxation oscillators to approximate the phase map near 
the singular limit and use this construction to predict the 








shape of the (finite) phase response curve where, for the 
sake of illustration, we stick to impulses, that is, u(-) = 
aS(-), and square pulses of finite duration, that is, u(-) = 

u [!+(•) - l+(- - A)]- 

A. Singular phase map and isochrons 

A first step towards the prediction of singular (finite) phase 
response curves is the construction of the (asymptotic) phase 
map and isochrons for the system ([T]) in the singular limit. 
Their construction relies on the fast-slow dynamics ([^-Q 
and is illustrated in Figure]^ 

Since the singular periodic orbit 7° is a one-dimensional 
piece-wise smooth curve in K^, it is naturally parameterized 
in terms of a single scalar phase on the unit circle S^. As in 
the nonsingular case, the phase map will be chosen such that 
the phase variable linearly increases with time. 

We choose to associate the zero-phase reference position on 
the singular periodic orbit with the lower fold {x-, Z-), that 
is 0°(a:_, Z-) =: = 0. As jumps are instantaneous in the 

singular limit, all points of the (weakly) unstable critical fiber 
joining (x_,z_) to (5-|_(z_), z_) are also associated with a 
phase equal to zero. 

The phase 6 associated with a point {x, z) is the “normal¬ 
ized” fraction of (slow) time At needed to reach this point 
along the reduced dynamics Q flow from the reference initial 
condition. For a point (a:i, zi) on the upper branch, the phase 
will be given by 

0°(a;i,zi) := uj ° Ati. 

For a point {x2, Z2) on the lower branch, the phase will be 
given by 

0°(a;2, Z2) := UJ° At+ -f a;° At2 

where the first term corresponds to the flowing time on the 
upper branch (up to the upper fold) and the second term 
corresponds to the remaining flowing time on the lower 
branch. To simplify notation, it is convenient to denote by 
0 °(a:+,z+) =: 0 + = Lo^ At+ the phase associated with the 
upper fold (and all points of the (weakly) unstable critical fiber 
joining {x+,z+) to {b-{z+), z+)). 

The notion of singular phase map can be extended to any 
point {x, z) in the basin of attraction of the singular periodic 
orbit. Because, in the singular limit, any singular trajectory 
starting from [x, z) instantaneously jumps from its initial 
condition to a branch of the critical manifold, all points on 
the same vertical line (that is, with the same value of slow 
variable z) jumping to the same branch are associated with 
the same phase. 

• All points on the line z = (resp. z = z+) are 
associated with the phase 0 _ (resp. 0 +). 

• For points with a slow variable in the bistable range, the 
asymptotic phase 0 i of a point (xi, Zi) belonging to the 
basin of attraction of the upper (resp. lower) branch is 
thus given by the phase 0i of the point at the intersection 
between the line z = zi and the upper (resp. lower) 
branch of the singular periodic orbit 7°. 


0®(a'i, Zi) = + CJg Ati 



Fig. 3. Geometric construction of (asymptotic) singular phase map. The 
phase map associates with each point on the periodic orbit a phase which 
corresponds to the normalized time cj® At required to reach this point 
from the reference position (x_, z-). For points on the lower branch, it is 
convenient to measure the normalized time from 2 +) and to add the 
phase 0+ Ar-f. Because all points on a same vertical ray (in the 

bistable region) and converging to the same branch instantaneously jump on 
the branch in the singular limit, the asymptotic phase map associates them 
with the same asymptotic phase. In addition, other vertical lines (outside 
the bistable region) are associated with the same phase because these points 
converge in the same Ar (mod T^) to (x+, 2 +). 

• In addition, all points outside the bistable range that 
converge to the upper fold in the same time interval 
Ar (mod T°) as {xi,zi) are also associated with the 
asymptotic phase 6i. 

An elegant way to summarize the definition of the singular 
(asymptotic) phase map is 

9-+ z,0) (mod 27r), if (Cl), 

0++ a;° ^_(z+, z, 0) (mod 27r), if (C2), 

with 

(C1)=(x,z)gS(5“)UT+, 

(C2) = ix,z) G S(5“) UT_. 

where ip,{zo, Zr, 0) (with • standing for + or —) are functions 
that measure the time needed to travel along the critical 
manifold from the initial condition zq to flnial condition Zr 
and where B{S^) is the set of points that jumps to the stable 
branch 5“ of the critical manifold. 

Singular isochrons are thus vertical lines for values of z 
outside the bistable range and vertical rays for values of z in¬ 
side the bistable range. In the bistable region, vertical rays are 
separated by the repulsive branch of the critical manifold. 
The vertical ray and the vertical lines associated with the same 
phase join at infinity (see Figure]^. 

For constant inputs u(-) = u, the function tp, (zq, Zt, u) can 
easily be computed by integrating the reduced dynamics Q 
on the stable branches of the critical manifold and they read 

1 

'lp,{zo,Zr,u)= ——- 
















Fig. 4. Effect of positive impulses in the fast-slow dynamics 0-0. 
(Case 1) Close enough to the lower fold (on the lower branch), the reset 
state crosses the separatrix (red curve) and converges toward the upper 
branch instantaneously. The phase shift corresponds to the phase difference 
corresponding to the skipped portions of the singular periodic orbit (green). 
(Case 2) Fai‘ from the lower fold (on the lower branch) or on the upper 
branch, the reset state converges back to the initial state instantaneously. As 
a consequence, no phase shift is produced. 


Remark 2: For presentation convenience, we intentionally 
do not consider the unstable branch of the critical manifold 5’' 
as being part of the basin of attraction of the singular periodic 
orbit. For small e, this repulsive branch is perturbed into a 
repulsive set which has zero Lebesgue measure. 

Remark 3: For convenience, the singular periodic orbit 7 ° 
is parameterized by the map x''' : —?► 7 ° that associates 

with each phase 6 on the unit circle a point (x'^(0), z'^(d)) on 
the periodic orbit. 

B. Singular (finite) phase response curves 

We derive the singular (finite) phase response curve for two 
inputs: impulses, that is, u(-) = a6{-), and square pulses of 
finite duration, that is, u{-) = u [!+(•) — !+(• — A)]. 

I) Impulse: An impulse u(-) = a5(-) induces a jump 
of the fast variable x in the fast-slow dynamics ([^-Q. The 
singular (finite) phase response curve is thus given by 

Q°(0; a<)(•)) = 0°(a;^(0), z'^{0) + a) - 0. 

As illustrated on Figure if the impulse lets the state 
cross the unstable branch of the critical manifold (case 1 ), it 
produces a phase shift. In the opposite case (case 2), the state 
converges back to the initial condition almost instantaneously. 

For simplicity, we assume monotonicity of this separatrix 
in the bistable region (that is, {dbr/dz){z) > 0). 

Given a positive impulse of amplitude a, there exists a 
critical value Zc{oi) of the slow variable such that a trajectory 
starting on the lower branch crosses the separatrix under the 
effect of the impulse for all z, such that z- < z < Zc{a). The 
critical value Zc{a) is given by 

Zc{a) = {z S K : ^-(z) + a = br{z)}. 

The asymptotic phase associated with this critical point 
{b-(Zc{o!)), Zc{a)) on the stable branch is denoted 6c{a). 


The phase shift A9 induced by an impulse corresponds to the 
portion of singular periodic orbit skipped due to the impulse. 
The phase response curve is given by 


Q°{e-ad{-)) 


9.+u:°'iP+{z_,z'r(e),O)-0, if(C3), 
0 , o/w, 


where (C3) stands for 9c{a) < 9 < 9-. 

Following a symmetric reasoning for negative impulses, 
that is, u(-) = —a6{-), the phase response curve is given by 




0++wO^_(z+,z^(0),O)-0, if(C4), 
0 , o/w, 


where (C4) stands for 9c{a) < 0 < 0+ and Zc(q:) = {z G K : 
6 +(z) -a = br(z)}. 

2) Square pulse of finite duration: A square pulse of finite 
duration uf ) = u [!+(•) — !+(• — A)] induces a behavior in 
the fast-slow dynamics @-(0 that is less trivial. 

The phase response curve is given by 


Q°(0; 7r(.)) = 0°(a:A(0), ^a(0)) - (0 + A°) (5) 


where (a:A( 0 ), -za( 0 )) is the state at time A^ for the reduced 
dynamics starting from (x^(9), z'’'(0)) where A° is the pulse 
duration in the slow time scale and in the singular limit. It is 
thus necessary to compute the state (xa,za) of the trajectory 
at the end of the pulse in order to compute the reset phase 
associated with its initial condition. 

In the following, we describe how we can compute the state 
(xa,za) using only the information contained in the func¬ 
tions 'tp-(z+ + u, z, u) and (see Figure]^. 

Starting from the initial condition (xo,zo) on the critical 
manifold, the trajectory evolves as follows (see Figure]^. 

(1) Under a constant input u, the critical manifold of the 
system is shifted along the z-axis. The singular trajec¬ 
tory jumps thus instantaneously to the branch of the 
“shifted critical manifold” corresponding to the basin 
of attraction to which the initial state belongs. 

(2) Then, the trajectory evolves on the “shifted critical 
manifold”, sliding slowly on branches and jumping 
instantaneously when it reaches “shifted folds”. 

(3) Finally, the trajectory jumps instantaneously back to the 
critical manifold at the end of the pulse. 

Because the slow variable z is one-dimensional, the evolution 
of a trajectory under constant input u on an attractive branches 
is fully characterized by the functions ' 0 -(z+ + u, z, u) and 
tjj_^.(z- + u, z, u) during the flowing time. The total flowing 
time has to be equal to the duration Aj. 

In Figure]^ we differentiate between two cases. In case 1, 
the initial condition on the lower branch of the critical mani¬ 
fold jumps directly to the upper branch of the shifted critical 
manifold. In case 2, the initial condition on the lower branch 
of the critical manifold jumps on the lower branch of the 
“shifted critical manifold”. Case 1 produces larger phase shift 
than case 2 . 

Remark 4: The duration A is expressed in the fast time 
scale, that is, A| = A. In the slow time scale, the duration 
is given by A^ = eAj. We assume the duration of the 









case 1 


case 2 



Fig. 5. Effect of positive square pulses of finite duration in the fast-slow dynamics 0-0. The state of the trajectory starting from initial 

condition (a^o^ ^o) (under a pulse of duration A) is graphically determined using functions ip» in order to predict the phase response associated with this 
pulse. The effect of a positive pulse is to shift temporally the critical manifold along the 2 :-axis to the right. The singular trajectory stalling from (xq, zq) 
evolves as follows: (1) jumps instantaneously on the shifted critical manifold, then (2) evolves around the shifted hysteresis (for a duration A = Aa + Ac), 
and finally (3) jumps back to the initial critical manifold. The main difference between case 1 and case 2 is that during step (1) the trajectory converges 
to the opposite branch (with respect to the initial point) of the shifted critical manifold (in case 1) or to the same branch (with respect to the initial point) 
of the shifted critical manifold (in case 2). 


pulse Aj (in the slow time scale) do not tend to zero in the 
singular limit and thus that the duration Aj tends to inhnity. 
This assumption is motivated by the fact that the duration of 
the pulse is often a fraction of the period. So we may have 
limf_>o A| = +00 and lime_>.o Tf = +oo, and a finite ratio 
lime_>.o Af/T/ = C (with C 0 and C ^ oo). 

V. Application to a neural oscillator model 

We illustrate our geometric approach on a simple neural os¬ 
cillator model developed by FitzHugh [12] and Nagumo [13]. 
This model is a popular two-dimensional simplihcation of the 
Hodgkin-Huxley model of spike generation 

V = V — /3 — w + u 

T w = a — bw + V 

where v is the voltage variable, w is the recovery variable, and 
e := l/r is a small parameter. 

A. Phase response curves for impulses 

Figure illustrates the (hnite) phase response curve of 
the FitzHugh-Nagumo model for excitatory impulses u(-) = 
a6{-), with a > 0. The solid line is the geometric prediction 
computed in the singular limit. Dots represent the phase 
response computed through numerical simulations of trajec¬ 
tories of the model for different values of the parameter e. 

The singular phase response curve is equal to zero except 
in one region of the periodic orbit which corresponds to 
the region right before the initiation of the upper part of 
the periodic orbit for an excitatory impulse. In this region, 
an impulse advances the initiation of the upper part of the 



Fig. 6. Phase response curves for excitatory impulses: singular geometric 
prediction (solid line) and numerical simulations (dots). (Parameter values: 
a = 0.7, b = 0.8, 7 = 1, |q| = 1.5.) 

periodic orbit. The phase advance decreases monotonically to 
zero until the phase corresponding to the lower fold. 

For small values of e, the geometric prediction matches 
very well the numerical phase response curves. For larger 
values of e, the prediction still matches (qualitatively) the 
larger phase shifts arising before the lower fold but do not 
capture the small phase shifts arising before the upper fold. 

B. Phase response curves for square pulses of finite duration 

Figure |7] illustrates the (hnite) phase response curve of the 
FitzHugh-Nagumo model for excitatory square pulses of hnite 
























Fig. 7. Phase response curves for excitatory pulses of finite duration: 
singular geometric prediction (solid line) and numerical simulations (dots). 
(Parameter values: a = 0.7, b = 0.8, / = 1, |ti| = 0.25, A = O.IT.) 

duration. The solid line is the geometric prediction computed 
in the singular limit. Dots represent the phase response com¬ 
puted through numerical simulations of trajectories of the 
model for different values of the parameter e. 

The singular phase response curve is equal to zero except 
in two regions of the periodic orbit. The first region which 
exhibits the highest phase shifts corresponds to same region 
as for the impulse case. The phase shifts in this region fol¬ 
low a piecewise law; the breaking point in the phase shifts 
corresponds to the separation between initial conditions that 
continue to evolve on the shifted initial branch and those that 
directly jump to the opposite branch. The second region cor¬ 
responds to point close to the other fold (see case 1 and case 2 
in Figure]^. An excitatory pulse may delay the termination of 
the upper part. 

Once again, for small values of e, the geometric prediction 
matches very well the actual phase response curves. For larger 
values of e, the prediction matches qualitatively both non-zero 
regions of the phase response curve. 

The main difference between the phase response curve for 
an impulse and for a pulse is that a positive pulse may delay 
the termination of the behavior on the upper branch, while a 
positive impulse may not. 

VI. Conclusion 

In this paper we overcome the limitation of the infinitesimal 
phase response curve approach to singularly perturbed oscil¬ 
lators by studying geometrically the asymptotic phase map 
and the finite phase response curve in the singular limit. By 
persistence results, this analysis provide semi-analytic predic¬ 
tion on the qualitative isochrons structure and on the (finite) 
phase response curve shape for arbitrary inputs. This result 
is illustrated on a relaxation oscillator forced by impulses 
and pulses of finite duration, confirming the goodness of the 
approach. 

Future work will aim at extending this analysis to more 
complex singularly perturbed oscillators, like bursters [10], 


and to oscillator synchronization studies, linking this result to 
fast threshold modulation phenomenon [14]. 
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